function [psd,f,ENBW] = mypsd(x,fs,Navg,order,wind,overlap)
% [f,psd,ENBW] = mypsd(x,fs,Navg,order,window,overlap)
% 
% Computes Power Spectral Density using Welch's method of 
% overlapped segmented averaging of modified periodograms.
% Based on techniques described in "Spectrum and spectral density
% estimation by the Discrete Fourier Transform", G. Heinzel, et. al.
% 
% INPUTS:
%   x = time series input [amplitude]
%   fs = sampling frequency [Hz] <default 1.0>
%   Navg = desired # of averages <default 1> 
%      NOTE: bin spacing = lowest frequency = fs/N
%   window = type of window ('flat', 'hanning', etc.)
%      Suported Windows:
%         - none (rectangular, boxcar)
%         - Bartlett (triangle)
%         - Hanning <default>
%         - Hamming
%         - Blackmann-Harris (92dB)
%         - Kaiser (kaiser3, kaiser4, kaiser5)
%   overlap = fractional overlap between sections 
%      NOTE: use overlap = 'auto' to use reccommended 
%            overlap for chosen window. <default 'auto'>
%
% OUTPUTS:
%   f = vector of Fourier frequencies
%   psd = power spectral density [amplitude^2/Hz]
%   ENBW = equivalent noise bandwidth [Hz]
%      NOTE: Other useful spectra can be determined as follows
%         - Power Spectrum, ps= psd * ENBW [amplitude^2]
%         - Linear Spectral Density, lsd = sqrt(psd) [amplitude/sqrt(Hz)]
%         - Linear Spectrum, ls = sqrt(ps) [amplitude]
%
% Ira Thorpe
% Updated 2-14-10

%% parse input arguments
error(nargchk(1, 6, nargin, 'struct')) 
if nargin < 6, overlap = 'auto'; end
if nargin < 5, wind = 'hanning'; end
if nargin < 4, order = 0; end
if nargin < 3, Navg = 1; end
if nargin < 2, fs = 1.0; end

%% Detrend data
n = (0:length(x)-1);
fit = polyfit(n,x,order);
x = x-polyval(fit,n);
%% Set up FFT Kernel

% get ROV (needed to get kernel length)
switch lower(wind)
    case {'rectangular','flat','boxcar','none'}
        ROV = 0;    
    case {'bartlett','triangle'}
        ROV = 0.5;
    case 'hanning'
        ROV = 0.5;
    case 'hamming'
        ROV = 0.5; 
    case {'blackman','blackman-harris'}
        ROV = 0.661;
    case 'hft116d'
        ROV = .782;
    case 'kaiser3'
        ROV = 0.619;
    case 'kaiser4'
        ROV = 0.67;
    case 'kaiser5'
        ROV = 0.705;
    otherwise
        error('Unknown window fucntion');        
end


% determine kernel length
N = 2*min(floor([(0.5*length(x)/(Navg*(1-ROV))) 0.5*length(x)]));

%disp(' ');

j = 0:N-1;                % points for time series
f_res = fs/N;             % frequency resolution
t = (1/fs).*j;            % time values (used in linear regression)
f = (0:N/2).*f_res;       % vector of Fourier frequencies

% Definition of windows (from G. Heinzel, et. al.)
switch lower(wind)
    case {'rectangular','flat','boxcar','none'}
        w = ones(1,N);
        %NENBW = 1;
        disp('No Windowing');
    case {'bartlett','triangle'}
        z = 2.*j/N;
        w = zeros(size(z));
        for i = 1:length(z)
            if z(i) <= 1
                w(i) = z(i);
            else
                w(i) = 2-z(i);
            end
        end
        %NENBW = 4/3;
        disp('Window is Bartlett');
    case 'hanning'
        z = 2*pi/N.*j;
        w = (cos((z-pi)/2).^2);
        %NENBW = 1.5;
        %disp('Window is Hanning');
    case 'hamming'
        z = 2*pi/N.*j;
        w = 0.54-0.46.*cos(z);
        %NENBW = 1.3628;
        disp('Window is Hamming');
    case {'blackman','blackman-harris'}
        % minimum 4-term B-H window
        % or 92dB B-H window
        z = 2*pi/N.*j;
        w = 0.35875-0.48829.*cos(z)+0.14128.*cos(2.*z)-0.01168.*cos(3.*z);
        %NENBW = 2.0044;
        disp('Window is Blackman-Harris');
    case 'hft116d'
        z = 2*pi/N.*j;
        w = 1-1.9575375*cos(z)+1.4780705*cos(2*z)-0.6367431*cos(3*z)+0.1228389*cos(4*z)-0.0066288*cos(5*z);
        %NENBW = 4.2186;
        disp('Window is HFT116D');
    case 'kaiser3'
        z = 2/N.*j-1;
        alpha = 3;
        w = besseli(0,pi*alpha.*sqrt(1-z.^2))/besseli(0,pi*alpha);
        %NENBW = 1.7952;
        disp('Window is Kaiser, alpha = 3');
    case 'kaiser4'
        z = 2/N.*j-1;
        alpha = 4;
        w = besseli(0,pi*alpha.*sqrt(1-z.^2))/besseli(0,pi*alpha);
        %NENBW = 2.0533;
        disp('Window is Kaiser, alpha = 4');
    case 'kaiser5'
        z = 2/N.*j-1;
        alpha = 5;
        w = besseli(0,pi*alpha.*sqrt(1-z.^2))/besseli(0,pi*alpha);
        %NENBW = 2.2830;
        disp('Window is Kaiser, alpha = 5');
    otherwise
        error('Unknown window fucntion');
        
end

% compute window sums and ENBW
S1 = sum(w);
S2 = sum(w.^2);
NENBW = N*S2/(S1^2);
ENBW = NENBW*f_res; 
 
% Determine overlap
if ischar(overlap)
    if strcmpi(overlap,'auto'), overlap = ROV; end
elseif isfloat(overlap)
    if overlap > 1 || overlap < 0, error('Overlap must be between 0 and 1'); end
else
    error('overlap must be either a numeric value or "auto" for automatic');
end
%disp(['Overlap is ', num2str(overlap*100), '%']);

% Determine PSD 
%disp('Computing PSD...');
%tic 

Navg = floor(length(x)/N*(1+overlap));      % # of segments
%disp(['Number of averages = ', num2str(Navg)]);
%p = zeros(Navg,N);

for i = 1:Navg;   
    j_low = floor(N*(i-1)*(1-overlap)+1);                % lower index for segment i
    j_high = j_low+N-1;                                  % upper index for segment i
    seg = x(j_low:j_high);                               % segment i
    fit = polyfit(t,seg,1);                              % linear fit to segment i
    res = seg - polyval(fit,t);                          % residual from linear fit for segment i
    win_res = res.*w;                                    % windowed residual
    y = fft(win_res);                                    % DFT of windowed residual
    p(i,:) = 2*y(1:N/2+1).*conj(y(1:N/2+1))/(fs*S2);     % PSD for segment i

end
psd = mean(p,1);    % Average over all segments to get PSD for entire data set
%disp(['Finished. Elapsed time = ', num2str(toc),' seconds.']);
